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Abstract. We compare the classical (mean-field) dynamics with the quantum 
dynamics of atomic Bose-Einstein condensates in double- we 11 potentials. The quantum 
dynamics are computed using a simple scheme based upon the Raman-Nath equations. 
Two different methods for exciting a non-equilbrium state are considered: an 
asymmetry between the wells which is suddenly removed, and a periodic time 
oscillating asymmetry. The first method generates wave packets that lead to collapses 
and revivals of the expectation values of the macroscopic variables, and we calculate 
the time scale for these revivals. The second method permits the excitation of a single 
energy eigenstate of the many-particle system, including Schrodinger cat states. We 
also discuss a band theory interpretation of the energy level structure of an asymmetric 
double-well, thereby identifying analogies to Bloch oscillations and Bragg resonances. 
Both the Bloch and Bragg dynamics are purely quantum and are not contained in the 
mean- field treatment. 



1. Josephson Hamiltonian 

The Josephson effect is a paradigm of macroscopic quantum mechanics that first arose in 
the context of superconductors. Josephson P predicted that a coherent current / oc sin 
would tunnel between two superconductors separated by a thin layer of insulator if there 
was a difference in the macroscopic quantum phase between the order parameters in 
the two superconducting regions (see [2] for a review). Josephson-type effects have also 
been realized in superfiuid ^He [3] , superfluid ^He |1] , and most recently in Bose-Einstein 
condensates (BECs) formed in atomic vapours [3 El d [H] . Trapped atomic BECs are 
well suited to fundamental studies of macroscopic quantum mechanics because almost 
everything about them can be controlled to a very high degree e.g. shape of the trapping 
potential, interatomic interaction strength, type of measurements performed etc., see 
[H] for a general review. This means that a wide range of parameter regimes can be 
achieved in a single experimental setup. In this paper we are interested in comparing 
and contrasting the 'classical' regime where mean-field theory provides an excellent 
description and a more quantum regime where quantum fluctuations play a role. 

The ac Josephson effect, which is driven by a difference in chemical potential 
between two sides of a tunnelling barrier, can be realized in an atomic BEC by 
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trapping the atoms in a double-well potential. Numerous theoretical studies, e.g. 
[Tot [TT| [T2| [T3| [HI [151 HSl HZj, have analyzed this setup. Unlike the situation found in 
strongly interacting systems, such as superconductors or quantum liquids, the Josephson 
equations for BECs can be derived using the mean-field Gross-Pitaevskii theory for 
which the underlying microscopic hamiltonian is well understood. Labelling the two 
wells by I (left) and r (right), the Josephson equations governing the evolution of the 
macroscopic phase difference = (pr—'Pi and the atom number difference k = {Ni—Nr)/2 
between two weakly coupled BECs in a symmetric double-well potential can be written 

m 

■ Ec,^ , 4k /N' Ej 

(p = ^-k -\ , = — cos (1) 

^ ^l-AkyN^ ^ 

k = - ^^l-AkyN^siiKp (2) 

where the dot represents a time derivative and is the total number of atoms 
N = Ni + Nr- Following the notation used for superconductors, the parameters Ej and 
Ec are called the tunnelling and charging energies, respectively. The tunnelling energy 
determines the maximum current Ij = Ej/h. By conservation of the total number of 
atoms, the current obeys I = Nr = —Ni = Ij\J^ — 4:k'^/N'^ sin 0. The charging energy 
arises from interatomic interactions and, providing A; <^ A^, it can be evaluated from the 
chemical potential at static equilibrium as 



E, 



c 



dN, 



(3) 

Ni=N/2 

Expressions for Ej and E^ in terms of microscopic quantities can be obtained from the 
Gross-Pitaevskii theory in both the tight-binding and Thomas- Fermi regimes |18j . 
By identifying and k as canonically conjugate variables, Hamilton's relations 
: IdH . IdH 

imply that the effective hamiltonian governing the dynamics of the macroscopic variables 
takes the form [131 [H] 

H ='^e - Ej^l- 4fc7Ar2 cos 0. (5) 

In this paper we will concentrate on the regime where the atom number difference 
between the wells is always much smaller than the total atom number /c -C iV. 
Furthermore, we assume the parameters obey E^ ^ Ej/N"^ which can always be satisfied 
for a large enough total atom number providing Ec^ Under these circumstances the 
hamiltonian (|5| reduces to [13] 

Hj = ^e-Ejcos<P . (6) 

and this is the form we shall work with from now on. This hamiltonian is analogous to 
that of a pendulum or, equivalently, to that of a classical particle moving on a sinusoidal 
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"washboard" potential. In the (former) latter case is the (angular) displacement and 
k the (angular) velocity. 

The relative phase and number difference k that appear in ^ are classical 
variables in the sense that their values are simultaneously well defined. This is in 
accordance with the derivation [13] of the Josephson equations ([T]) and ^ from the 
Gross-Pitaevskii theory which is a mean-field theory that assumes all atoms share 
the same macroscopic wave function (the condensate order parameter). The Gross- 
Pitaevskii theory has proved enormously successful as a description for atomic BECs 
trapped in single-well potentials. However, the double-well system provides a very 
simple and analytically tractable extension in which we can explore beyond mean- 
field effects essentially because the single-particle kinetic tunnelling/hopping energy 
represented by Ej can easily be much smaller than the interaction energy E^. Under 
these circumstances it is necessary to (second) quantize the Josephson hamiltonian ([6]). 
We do this by promoting and k to operators which satisfy [9] 

[0, k] = i. (7) 

In the 0-representation where k = — id/d0 the quantum version of the Josephson 
hamiltonian is 

E 

Hj = - Ejcos(f) . (8) 

To complete the quantization we also need to stipulate that the wave function -0(0) that 
the hamiltonian ^ acts on is single valued, i.e. is periodic -0(0 + 27r) = -0(0). This 
ensures that the eigenvalues of k are integers and is in accordance with the notion that 
is a phase. The difference between the dynamics generated by the classical ^ and 
quantum (|8| hamiltonians is the main theme of this paper. 

The requirement that the wave function be 27r-periodic is very natural from the 
point of view of the pendulum analogy. However, from the point of view of the 
particle in a washboard potential analogy it is a more restrictive condition. The time- 
independent Schrodinger equation associated with the hamiltonian (|8| is the Mathieu 
equation {Hj — E)"^ = |19j. According to the Floquet-Bloch theorem the general 
solutions of the Mathieu equation can be written -0(0) = exp(ig0)[/g(0) where q is the 
quasi-momentum. However, because our wave function is 27r-periodic we must set the 
quasi-momentum q to zero. Thus, it appears that analogues of a number of phenomena 
familiar from the physics of waves in periodic potentials, such as Bragg scattering and 
Bloch oscillations, must be absent from the Josephson problem because they require 
finite values of the quasi-momentum. On the contrary, we shall see in Sections [7] and 
[8] that there is a sense in which we can achieve finite q values and hence realize Bragg 
scattering and Bloch oscillation analogues in double-well systems. 

The validity of the quantization procedure given above to obtain the hamiltonian ^ 
is actually far from obvious |20]. The original iV-atom double- well system corresponds to 
a quantum many-body system which is then approximated by a Gross-Pitaevskii mean- 
field theory to give the Josephson equations ([T]) and The system is then re-second 
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quantized by quantizing the mean-field theory to give the hamiltonian ([8]). However, it 
turns out that in the regime k <^ N and Ec ^ Ej/N'^ the hamihonian ([s]) agrees with 
that obtained from a treatment based upon the fully quantum Bose-Hubbard model, 
see, e.g. [2T| . 

The organization of the rest of this paper is as follows. In section |2] we introduce 
the Raman-Nath equations which provide a simple framework for calculations associated 
with the quantum Josephson problem. In sect ions |3] and |4] we generalize our treatment to 
include an asymmetry between the two wells. In section |5] we compare the quantum and 
classical dynamics of the Josephson junction following an excitation created by taking a 
system which is at equilibrium in an asymmetric double-well and suddenly removing the 
asymmetry, which is the current standard experimental probe. Section [6] treats double- 
wells that have an asymmetry that is modulated periodically in time and Section [7] 
examines tunnelling resonances that occur at certain values of the asymmetry that are 
analogous to Bragg scattering. In Section [8] we present an analysis of the asymmetric 
double-well problem based on band structure theory and relate adiabatic sweeps of the 
asymmetry to Bloch oscillations. We also suggest a way to generate Schrodinger cat 
states. 



2. Raman-Nath equation for a BEC in a double-well potential 

The 27r-periodic wave function t) which determines the values of the macroscopic 
variables and k can be expanded as a Fourier series 

oo 

*(0,t) = -= J2 ^n(t)exp(in0) (9) 



where the 'plane wave' basis states exp (in0) / v 27r are eigenfunctions of the number 
difference operator k with integer eigenvalues n. The probability amplitudes An obey the 
normalization J2n l^nP = 1- In the original problem with N atoms (where without loss 
of generality we take to be even) the integer n must lie in the range —N/2 < n < N/2, 
but since we are working in the regime where the number difference is always small in 
comparison to the amplitudes An become vanishingly small long before \n\ = N/2 
and so we have extended the sum in ^ to ±oo. Substituting ^ into the Schrodinger 
equation associated with the Josephson Hamiltonian ([sj) 

Ec 



ih^^^{<f),t) 



- Ej cos(0) 



vl/(0,t) (10) 



2 ( 

yields the coupled infinite set of Raman-Nath (RN) differential-difference equations 

i^Anir) = n'^Anir) - ^ [A„+i(r) + A„_i(r)] . (11) 
clr Z 

We have written the RN equations in dimensionless form by defining r = Ect/2h and 

the ratio of the energy parameters as A 

A-^. (12) 
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The RN equations first arose in the context of the dynamical diffraction of light by 
ultrasonic waves in fluids |22l [23] and have subsequently found extensive use in the 
description of the diffraction of atoms by standing- waves of light [2lj. By numerically 
integrating in time a suitably truncated set of RN equations they provide a simple 
scheme for computing the dynamics of the macroscopic variables of the double-well 
problem. A novel feature of the atomic BEC realization of the Josephson junction 
is that great control can be exerted over A. Either by adjusting the intensity of the 
laser which forms the central tunnelling barrier, or by using a Feshbach resonance to 
manipulate the interactions, the magnitude of A can be varied between essentially zero 
and infinity. These experimental 'knobs' can also be turned during the course of an 
experiment leading to a time-dependent A which can also be easily handled within the 



RN framework (11). 



An alternative method for describing the dynamics of the macroscopic variables 
of the double- well problem is to expand \E'(0, t) in terms of the eigenf unctions of the 
Josephson hamiltonian. Let 

vi>(0,t) = J2a,^^i^,t) = J2a,i;^i<P)exp{-ieH/h) (13) 
j j 

where ip-'{(f)) is the jth eigenfunction and has an energy . Defining the scaled energy 
= 2e^ /Ec, the eigenstates obey a Helmholtz equation which has the same general 
form as the Mathieu equation 

"d02 

Similarly to above, we can expand the eigenfunctions in a number state basis 



A cos(0) 



ipj = E^ip^. (14) 



n=oo 

^'■(0) = ^ E Aiexpiin^) (15) 



n= 



which leads to the coupled set of time-independent RN equations 

E^Ai = n'Ai - ^{Ai^, + AU). (16) 

The time- independent RN equations thus take the form of recurrence relations describing 
a tridiagonal matrix whose jth eigenvalue E^ corresponds to a column eigenvector 
{. . . ^^3, A-'_2, ^-1, Aq, A{, A2, A3 . . .} made up of specific values of the amplitudes A^. 
It is, of course, exactly the recurrence relation obeyed by the Fourier components of 
the even and odd Mathieu functions cer{(p) and 56^(0), respectively |T9]. We shall 
stick to the notation ip^^cj)) to cover both cer{4>) and 56^(0) so that '0°(0) = ceo(0), 
= sei(0), ■0^(0) = cei(0), -0^(0) = 562(0) etc. As is to be expected from a 
problem involving a one dimensional wave equation, the eigenfunctions alternate in 
parity as one goes up in energy (i.e. increases the index j), with the ground state even. 
The Fourier space eigenvectors {. . . A-'_^, Aig, ^-i, A-f^, A{, A2, A^ . . .} are localized in k- 
space, i.e. — >■ for large enough n. In fact, the decay of A{^ with n is exponentially 
fast for large n [25]. Comparison with the tridiagonal matrix representing the exact 
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Bose- Hubbard hamiltonian [21] shows that corrections to (16) are of order 1/A^ and so 
become neghgible as oo. 

From Equation ( 14 ) we see that in the scaled variables the parameter A determines 
the height of the sinusoidal potential and thus the separatrix for the classical motion. 
The separatrix divides classical phase space into two qualitatively different types of 
motion. When the energy is smaller than A then in the classical pendulum analogy we 
have librational motion meaning that the phase (f){t) and angular velocity k{t) oscillate 
and so periodically reverse their signs. In the classical particle in a sinusoidal potential 
analogy this corresponds to the particle having an energy less than the barrier tops 
and so rolling around within a single well of the sinusoidal potential. In the physical 
double-well problem librational motion corresponds to Josephson plasmon excitations 
[8j . When is greater than A then in the classical pendulum analogy we have rotation 
meaning that the phase continuously winds up in only one direction (the phase lives 
on a torus with = —it and (p = it identified so that it lies in the range — vr < < vr) and 
k does not reverse its sign. In the classical particle in a sinusoidal potential analogy the 
particle has enough energy to roll into the neighbouring wells of the sinusoidal potential. 
In the physical double-well problem rotational motion is known as macroscopic quantum 
self-trapping [13] and corresponds to motion in which the system is locked in a state 
with a larger number of particles in one well despite the symmetry of the potential. 
However, quantizing the system means that quantum tunnelling though the sinusoidal 
potential barriers [E^ < A), and quantum reflection above the barrier [E^ > A) blurs the 
distinction between libration and rotation. We thus expect motion near the separatrix 
to be one place where quantum effects are particularly visible. 

As A is the only parameter in our hamiltonian, its magnitude plays an important 
role. In the literature three regimes are usually identified: (1) the Rabi regime A ^ A^^, 
(2) the Josephson regime > A > 1, and (3) the Fock regime A ^ 1. In particular, 
the Rabi and Fock regimes correspond to the non-interacting and interaction-dominated 
limits, respectively, and the Josephson regime lies in between, see [TTl [9l |2T] for more 
discussion. In this paper our assumption E^ ^ Ej/N"^ excludes the Rabi regime. In 
fact, for the most part we shall be in the Josephson regime, with the exception of the 
discussion of the Bragg scattering and Bloch oscillation analogies in Sections [7] and 
[8] which concern the border between the Josephson and Fock regimes where A < 1. 
Having reduced the problem to just two regimes, from now on we take the view that 
A plays a role analogous to the dimensionless ratio of the classical action to Planck's 
constant. More precisely, A = 2Ej/ Ec = {S/hY, where is a constant having the units 
of action. The comes from the kinetic energy E^. With this identification we see that 
A determines how quantum the system is: 

• A small. In this case the system may be viewed as being in a very quantum regime 
in the sense that there are only a few quantum states below the separatrix, i.e. 
'trapped' in the sinusoidal well. Only a small tridiagonal matrix (16) is required to 
capture the states having energies up to the separatrix and the classical Josephson 
equations ([l| and ^ are expected to give a poor description of the dynamics below 
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the separatrix. 

A large. When A is large we enter the semiclassical limit where there are many 
eigenstates inside the sinusoidal well and a large tridiagonal matrix is required. 
The classical Josephson equations are expected to give a good description of the 
dynamics (except close to the separatrix). 



In the semiclassical limit a large tridiagonal matrix (16) is required to capture the 
eigenstates up to the separatrix. An estimate of the required minimum dimensions 
A/" X A/" of the matrix size is given by [25] 



J\f = V2A . (17) 

Eigenstates well below the separatrix (which are localised near the bottom of the 
sinusoidal wells) have a linear energy spectrum like the harmonic oscillator (see Figure 
[T|. These states are known as Josephson plasmon excitations and occur at integer 
multiples of the energy f9j 



hupi = ^E,Ej . (18) 

As shown in the inset in Figure [T]^a), above the separatrix the eigenstates rapidly tend to 
degenerate pairs. Indeed, well above the separatrix (as j oo) the sinusoidal potential 
becomes irrelevant and the hamiltonian tends to that of the quantum rotor. One of 
the eigenstates in each pair is ccj and has even parity and the other is scj and has odd 
parity: they are approximately the (±) superpositions of clockwise and anticlockwise 
rotor states. Nevertheless, the sinusoidal potential does lead to a small energy splitting 
between the two states of each pair that scales as E^^^ — E^ = 0{A^ / as j — >• oo 
P^ . Taking each pair as a single unit, the spectrum of the units is quadratic as j ^ oo 
as expected for the quantum rotor. 

The energy spectrum flattens out near the separatrix (see Figure Q. This can 
be understood in terms of the divergence of the period of the classical motion at the 
separatrix which in turn causes the density of states (oc 1/period) to have a peak 
there [20] ■ The density of states can be calculated from the numerical eigenvalues as 
D{E^) = 1/\E^~^^ — E^ and analytic expressions can be derived using Bohr-Sommerfeld 
quantization. The expressions valid below and above the separatrix are, respectively 

[SHEHIEE] 




where K(x) is the complete elliptic integral of the first kind [I9]. As can be seen in 
Figure [TJ the numerical and analytic expressions are in excellent agreement except right 
at the separatrix where the Bohr-Sommerfeld method breaks down. Although we shall 
not make use of them here, analytic solutions to the RN equations in the semiclassical 
limit are available [25] . These are based on uniform approximations that are valid right 
through the separatrix and so go beyond WKB/Bohr-Sommerfeld quantization. 
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Figure 1. A = 25000. (a) The energy eigenstates versus the eigenstate number j. The energy 
spectrum is Unear for energy states far below the classical separatrix (which lies inbetween eigenstates 
285 and 286). (b) The density of states as a function of energy. We show both the exact numerical 
result and the approximate analytic result given in ( |19| . The agreement is excellent except very near 
the separatrix. Both predict a sharp peak in the density of states at the separatrix. 



3. Asymmetric double-well potential 

In order to investigate Josephson oscillations in the double-well potential it is necessary 
to first excite the system into a non-equilibrium state so that its subsequent dynamics 
can be observed. One way to do this is to start from an equilibrium state in a slightly 
asymmetric (tilted) double-well potential and then to suddenly make the potential 
symmetric. The two key experiments [7] and [8] have both used this method and we 
shall model this situation in this section. 

The Gross-Pitaevskii equation for a BEG in an asymmetric double-well leads to the 
following Josephson equations [13] 



4k 1 

h(j) = Eck + Ej— , =cos0 + Ae, (20) 



hk = - Ej^l-4k'^/N^sm(f) (21) 
where Ae is the difference between the zero-point energies of the two wells, i.e. magnitude 



of the tilt. The equilibrium state is defined by = and = 0, and so from (21 ) we see 
that the equilibrium phase difference in the asymmetric potential is still zero: 0eq = 0. 
However, there will be an unequal number of atoms on the two sides, i.e. keq 7^ 0, with 
more atoms sitting in the lower well. When the potential is suddenly changed to being 
symmetric our initial conditions are therefore 0|t=o = and k\t=o 7^ 0. In the pendulum 
analogy this corresponds to the pendulum starting at the instant where it is pointing 
vertically downwards but with a finite angular velocity. 



The equations of motion (20) and (21) are generated by the classical hamiltonian 



En I Ak"^ 

Ha=^k^-EjJl-—cos<P + Aek . (22) 
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Figure 2. A = 25000. Here we plot the equilibrium value of the relative population imbalance as a 
function of the energy asymmetry. The exact result from solving | |24| l numerically is compared with the 
approximate result given by Equation j26| . (a) The number of particles is A'^ = 10'^. The agreement 
is reasonable until AE ^ N. (b) A'^ = 10°. The agreement is extremely good until AE = N. 



In the regime k/N <^ 1 and ^ Ej/N"^ this takes the simphfied form 

E 

= -^k^ - Ej COS0 + Ae k 



(23) 



which will be referred to as the asymmetric Josephson hamiltonian. If the asymmetry 
is too large we risk violating the condition k/N <^ 1 even for the equilibrium state. It 
is therefore important to establish this extra condition of validity upon the hamiltonian 



(23). At equilibrium Equation (20) becomes 

2 



^'^eq 



/AE 



m \ 2 



+ h 



eq 



4A^ 



( ^eq \ 



(24) 



where AE = 2Ae/Ec is the dimensionless tilt asymmetry parameter. Assuming that 
keq/N is small we can expand as 

1 



AE 



k 



eq 




4fce\ 



2 m 



+ 



-2A^ 

iV2 



This gives 



k 



eg 
N 



AE N 
'^2A + A^2 



(when kejN < 1) 



(25) 



(26) 



In Figure |2] we compare (26) with the exact result obtained by numerically solving (24). 
We see that the two are in excellent agreement almost all the way up to AE = N 
which is the saturation point where all atoms have moved into a single well. Our 
earlier assumption A^^ ^ A (exclusion of the Rabi regime) means that we can further 
approximate (26) as 



k 



AE 



eq 



(27) 



We therefore see that the extra condition that the pendulum hamiltonian (23) is valid 
in the asymmetric case is that AE -C A^. Note that the result (27) is actually the exact 
prediction given by the pendulum hamiltonian (23). 
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Figure 3. A = 25000. The ground state probability densities (as a function of population difference) 
for three different asymmetric double-well potentials: AE = 200,100, and for (a), (b), and (c) 
respectively. For zero asymmetry the expectation value of m is zero meaning the analogue quantum 
pendulum will remain motionless (apart from zero-point fluctuations). When the asymmetry takes 
some finite value the system has the greatest probability of being found in a state where m is equal 
to —AE/2, which is the same as the classical prediction ||27|l. 



4. Raman-Nath equation for an asymmetric double-well 

The wave function for the asymmetric double-well is still 27r-periodic and so for a given 
value of AE we expand the jth eigenstate of the system as 

^i = ^T.BLexpitm<P) (28) 
V2vr m 

where the subscript a denotes "asymmetric". The time-independent Schrodinger 



equation with the Hamiltonian (23) gives 



EiBi^ = im' + AE m)Bl - ^{bI^, + i?Li)- (29) 

Figure |3] plots the ground states for various values of AE. An asymmetry produces 
a non-zero expectation value for the population imbalance, as we saw classically. 
Suddenly switching off the asymmetry [2T] propels the system into motion and from 
the perspective of the macroscopic quantum mechanical variables we assume that this 
process can be modelled by a projection of the equilibrium quantum state in the 
asymmetric potential onto the eigenstates of the symmetric potential. For this purpose 
it is useful to relate the two sets of eigenstates via the matrix of coefficients Cmn 

= E (^mn^"" ' (30) 



Expanding the eigenstates in the number difference basis like in ([15j) and (28) we find, 

Cmn = ii^nra) = ^T. B^ f exp[i(g - p)m = E ApB; (31) 

where we have used the fact that the amplitudes are real (as are B^). For simplicity 
we take the initial state to be the ground state in the asymmetric potential. The resulting 
projection coefficients Cji are plotted in Figure |4] for different initial asymmetries AE. 
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Figure 4. A = 25000. The initial state expressed in terms of the eigenstates of the symmetric 
double-well hamiltonian for 3 different energy shifts, i.e. immediately following 'symmetrization' of 
the wells, (a) AE = 400 — > -E^, = 80.2. This corresponds to exciting a set of eigenstates entirely 
below the separatrix. (b) AE = 445 — > Ex = 99.2. This corresponds to exciting a set of eigenstates 
in both regimes centered around the separatrix. (c) AE = 600 — > Ex = 180.4. This corresponds to 
exciting a set of eigenstates above the separatrix. 



To help gauge the degree of excitation generated by each value of AE we introduce the 
notation E^. This is the expectation value of the excitation energy in the symmetric 
double- well expressed as a percentage of the separatrix energy. Thus E^ = corresponds 
to the ground state energy of the symmetric double-well and E^ = 100 to the separatrix 
energy. We see from Figure |4] that when exciting below the separatrix the distribution 
is smooth and gaussian-like, roughly corresponding to a coherent state. Excitations 
above the separatrix are no longer smooth but oscillate strongly. These two distinct 
behavioural regimes are joined at the separatrix which has properties of both. 



5. Classical versus quantum dynamics 



We now consider the dynamics of the macroscopic variables and k following excitation 
by the method described in Section |3} The classical (mean-field) dynamics are governed 
by Josephson 's equations 

d(f) 

&r 
dk 

d7 

with initial conditions 0(0) = and A;(0) 



2k. 



Asini 



(32) 
(33) 



-AE/2, see Equation (27). Analytical 



solutions to (32) and (33) with the specified boundary conditions can be found in terms 



of special functions. For example, when the motion is below the separatrix we have 



= - 2arcsin [sn i^AE r/2 | 8A/(AE)2 
k = - {AE/2) dn (AE r/2 | 8A/(AE)2 
where sn(^|m) and dn(6'|m) are Jacobian elliptic functions [19 



(34) 
(35) 

When the motion is 



above the separatrix Equation (35) for k{t) remains the same but Equation (34) for 0(t) 



must be adjusted so that when reaches either ±7r it is then mapped to =F7r so that 



the evolution on the phase torus is continuous. The solutions (34) and (35) are periodic 
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with period 



To 



{8/AE) g?{K (8A/(AE)2)} 



(36) 



where 3? denotes the real part. Above the separatrix the motion is still periodic but 
the expression (36) for the period must be divided by two. When A/(Ai?)^ ^ 1 we 
can asymptotically expand the complete elliptic integral K(m) to find the period of the 
low-lying Josephson plasmon excitations 



To ~ TT 




lAE^ 9 /AE^^ 
^ 4 8A ^ 64 \ 8A ^ 



+ 



when 



A 

AE^ 



oo 



(37) 



where the first term corresponds to the harmonic approximation. The opposite limit, 
namely A/AE"^ <^ 1, is relevant for the high- lying rotor excitations: well above the 
separatrix we have 



T"0 



27r 
AE 



1 + 



1 8A 



4AE2 



9 / 8A y 
^ 64 [aE^J 



when 



A 



AE^ 



< 1 



(3J 



The quantum dynamics are treated using the RN equations which we use to 
calculate the expectation values of the operators and k. Whether the time-dependent 
version (11) or the eigenfunction version (16) of the RN equations is more suitable 
depends upon the length of time we want to track the dynamics for. For short times it 
is more efficient to use (11), but for longer times (16) is in principle faster because in 
this case time evolution is accounted for purely by the phase factors attached to each 
eigenfunction. Starting with the eigenfunction version, we expand \l/(0, t) in terms of 
eigenf unctions as in (13) and find 

(0(^)) 



(^( 



, T 



v[/(0,r)) 
i{E^ 



Substituting in the Fourier series (|15|) for the eigenfunctions we obtain 



(0(r)) 



n — m 



'ajA'^Ai sin 



E E 

j,k m,nj^m 

A similar calculation for the expectation value {k{T)) yields 
{k{T)) = V m(a,.Ya,AtAL cos \(E'' - E^)t 



J2 m{akT a, A^Al cos (E' 

j,k,m 



(39) 



(40) 



(41) 



To compute (40) and (41) we need to know the coefficients aj of the eigenfunction 
expansion (13). These are precisely the coefficients (31), i.e. aj = Cji for the case where 
the initial state is the ground state in the asymmetric double-well. 

Turning to the time-dependent version of the RN equations as given by (11), the 
time- dependent amplitudes A„(r) are evolved from their values at r = which are 
directly given by those of the ground state in the asymmetric double-well: ^^(0) = Bj^, 
see Figure |3} Using ^ the expectation value of the relative phase is 

(0(r)) = (v[/(0,r)|0|^((/.,r)) 
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'^'^ n,m •'-'^ 

^ _ A:(r)A^(r)(-l)— " ^^g^ 
m — n 

The time-dependent amplitudes A„(r) are in general complex numbers, so that the 
above expression is always found to be real. The expectation value {k) is calculated 
similarly, 

(fc(r)) = -(vl/(</>,r)|i^|vl/(<^,r)) = j:^\A.^i^)\' . (44) 

We have verified that the two versions of the RN equations give identical predictions for 
the expectation values. However, note that the summations appearing in (43) and (44) 
are over fewer indices than the equivalent summations in (40) and (41). This tends to 
make computations based upon the time-dependent version of the RN equations faster, 
especially when working in the semiclassical limit where the number of amplitudes that 
need to be included becomes large. 

In Figures |5| [6} and [7] we compare the classical and quantum predictions for 
(f){t) and k{t) for three different regimes of excitation: far below, near, and far above 
the separatrix, respectively. All the upper graphs plot the temporal evolution of the 
population imbalance and the lower graphs plot the temporal evolution of the relative 
phase. Like in Section [4], the degree of initial excitation is specified by E^. For 
excitations that are far below the separatrix we see that the quantum and classical 
predictions agree well for large A even for times corresponding to many classical 
periods but when A becomes small, corresponding to a more quantum system, there 
are observable differences in the expected frequency of oscillation. This renormalization 
of the classical frequency by quantum fluctuations has been discussed previously by 
Smerzi and Raghavan [22] ■ In either case the classical period far below the separatrix 



is accurately given by the flrst few terms in the expansion (37). 

For excitations very near the separatrix the classical prediction is valid only for short 
times for any value of A. This is because the classical motion is qualitatively different 
above and below the separatrix whereas the quantum motion contains elements of both 
due to the combined effects of quantum tunnelling and the flnite width in energy of 
the initial wave packet (see Figure |4]). We see in Figure |6] that near the separatrix 
the quantum and classical predictions diverge from each other on a time scale that 
is always shorter than one quarter of a classical period tq (the quantum result clings 
longest to the classical one as A — > oo). This can be understood by noting that the 
initial state has (0) = and is localised around the bottom of the sinusoidal well 
where quantum and classical agree best. However, for motion below the separatrix 
(as in Figure |6]) the subsequent evolution always reaches the classical turning point 
0tp = — arccos[l — A£'^/(4A)] at times equal to one quarter of the classical period: 
Ttp = To/ 4:- At the classical turning point the classical motion reverses direction but 
part of the quantum wave packet tunnels through the barrier (which is narrow near 
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(a) A = 10 (b) A = 150 (c) A = 25000 



Figure 5. Comparison of the classical (dashed curves) and quantum (solid curves) predictions for 
the temporal evolution of the population imbalance (top row) and relative phase (bottom row) for 
three different values of A (which is related to Planck's constant as A oc h~^). These dynamics are 
for excitations far below the separatrix. (a) AE = 0.267 Ex = 0.1. (b) AE = 1 ^ Ex = 0.1. (c) 
AE = 10 ^ Ex = 0.1. 





2 4 6 1 2 0.04 0.08 



(a) A = 10 (b) A = 150 (c) A = 25000 



Figure 6. The same as Figure [s] but for excitations near the separatrix. (a) AE = 8.44 Ex = 
99.9. (b) AE = 34.141 ^ Ex = 99.9. (c) AE = 445 ^ Ex = 99.2. Even in the semiclassical limit 

A ^ oo the classical and quantum results only agree for the first 1/4 period. 



the separatrix) with the result that the quantum and classical predictions diverge at 
or slightly before this point. For the parameters chosen in Figure [6] neither of the 
expansions (37) and (38) provide particularly accurate approximations to the exact 
result (36) for the period tq but (38) gives the right order of magnitude indicating that 
rtp = 0{7c/AE). 

Excitations far above the separatrix are rotating states in the classical pendulum 
analogy. In Figure [7] we can clearly see the macroscopic self-trapping effect in the 
behaviour of the number difference k in both the quantum and classical predictions: 
the amplitude of k performs small (relative) oscillations about a particular fixed value 
and does not reverse sign. The striking difference between the quantum and classical 
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Figure 7. The same as Figure [s] but for excitations far above the separatrix. (a) AE = 18.9 — ► 
Ex = 501.1. (b) AE = 76.4 E^ = 501.1. (c) AE = 1000 ^ E^ = 501.1. Note the change of scale 
on the abscissa of the (c) panels: the oscillations become very rapid as a function of t in the classical 
limit and so we have only shown the first few. 



predictions is that the quantum prediction undergoes periodic collapses and revivals due 
to the discrete energy spectrum. In general, the time scale for revivals of the quantum 
wave function is given by [30] 



^- ^ WW\ 

where E"{jo) is the 2nd derivative of the energy spectrum with respect to the quantum 
number j labelling the energies and is evaluated at the centre of the wave packet jo- 
Far above the separatrix the quantum Josephson hamiltonian reduces to that of the 
quantum rotor 

which has the spectrum = Ecj'^/2 where j is an integer j = 0, ±1, ±2, . . .. Thus, 
in our scaled time units, we obtain Trev = 27r. This is actually twice the revival time 
of 7r that can be clearly observed in Figure [7](a). The discrepancy can be explained by 
noting that in Figure [7] we have plotted the expectation value which contains the square 
of the wave function and thus we expect the expectation value to revive on a time scale 
which is half that of the wave function revival time |31|. From Figure [T] we also see that 
the revival time increases as the system becomes more classical in the sense that more 
oscillations occur during the time Tj-ev We note in passing that collapses and revivals 
can in principle also take place for wave packets excited below the separatrix but the 
time scale is typically much longer than above the separatrix. This is because below 
the separatrix the spectrum is to the first approximation linear and for a purely linear 
spectrum then Equation (45) predicts T^v —>■ oo [30j. Thus, only the small non-linear 
correction terms contribute to a finite collapse and revival time for excitations far below 
the separatrix. 

As a final point in this section we check that the uncertainty relation a^^ak > 1/2 
satisfied by the variances cr^ and ak of the phase and number difference is obeyed for the 
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Figure 8. These 3 curves demonstrate the uncertainty relation for the dynamical variables <j> and 
k. (a) AE = 1 — > Ex = 0.1. For excitations far below the separatrix there are small fluctuations, 
but the uncertainty remains very near 1/2. (b) AE = 34.141 — > Ex = 99.9. For excitations near the 
separatrix, the uncertainty begins at a minimum and then increases sharply after a short time, (c) 
AE = 100 — > Ex = 857.9. For excitations far above the separatrix the uncertainty begins to oscillate 
but centers itself around a finite value. 



evolution we have computed. In Figure [8] we plot the uncertainty product for A = 150 
and for each excitation regime. Excitations far below the separatrix remain very close 
to the minimum uncertainty, experiencing only small fluctuations, as expected for a 
gaussian (coherent) wave packet in a nearly harmonic potential. Near the separatrix 
the uncertainty product starts at the minimum value of 1/2, then after a short time 
increases very rapidly. Far above the separatrix the uncertainty product starts as a 
minimum but then oscillates and centers itself around a finite value above 1/2. 

6. Time-modulated double-well potentials 

The method of excitation described in Section [3] generates 'classical-like' gaussian wave 
packets of eigenstates. However, to fully explore the quantum dynamics of the double- 
well system it is desirable to be able to excite individual eigenstates, for example to 
generate the Schrodinger cat states discussed in Section [8] below. One way to excite 
an individual eigenstate is to start from the ground state and apply a time-modulated 
asymmetry which is resonant with a particular transition. As an example we will use this 
method to excite a single eigenstate near the separatrix. States near the separatrix are 
interesting because they are at the boundary between two qualitatively different classical 
regimes. The special behaviour of these states can be hidden due to the population of 
the surrounding states if the method of Section [3] is used. 

To transfer the system from an eigenstate a of the symmetric double-well with 
energy E"- to another eigenstate b of the symmetric double-well with energy we 
apply a time-varying asymmetry of the form AE{t) = ({t) sin(i7r) which is resonant 
with the transition so that Q = {E^ — E°')/h. The pulse envelope C(''")) which is assumed 
to be slowly varying in comparison to h/Q, determines the temporal shape of the pulse. 
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Figure 9. Excitation of a single eigenstate (j = 5) from the ground state {j = 1) by means of a 
time oscillating potential. Upper panel: snap shots after various times of the probabilities \aj\'^ to 
occupy the jth eigensate of the balanced well system, (b) In order to reach the state j = 5, which lies 
immediately below the separatrix, we pump energy into the system in a series of four pulses (whose 
profile is of the form a X (1 + tanh[T/b])) and whose area is determined by | |48| l. Each pulse is resonant 
with its respective transition. We set the total time of the process to be r = 400. Notice that the 
final probability to be in state 5 is not exactly 1 because other states have been marginally excited. 
This result can be improved if we make the pulses more adiabatic, requiring smaller pulses over longer 
times but conserving the total area. 



The evolution of the amphtudes is governed by the time- dependent RN equation 

i^-J^ = in' + AE(r)n)A„(r) - -(A„+i(r) + A„_i(r)). (47) 
dr 2 

In what follows we shall choose the initial state to be the ground state of the symmetric 

double-well. 

A good estimate of the time required to excite the target eigenstate b can be 
obtained by approximating our multi-level system as a two-level system consisting 
of just the states a and b and considering the effect of applying a time-modulated 
asymmetry with a constant amplitude (. The characteristic frequency at which the 
two-level system is driven between its two levels is then given by the Rabi frequency 
a;^ = —i(Ec/ {2h)\{'ilj^\k\'ilj"')\ |32j. Introducing the dimensionless Rabi frequency 

= 2huj^^/E^ we find 

= CT.rnAlA':^ . (48) 

m 

The time to fully transfer the system from a to 6, i.e. to deliver a vr-pulse, is given by 

In practice, the amplitude C(r) must be switched on smoothly from zero so that 
the frequency spread is small and only a single upper eigenstate is excited. Empirically 
we find that a single target state can be excited only if the maximum value of C(r) 
obeys the conditions Cmax ^ A and Cmax -C Vt^. Nevertheless, by direct integration of 
the RN equations (47) using the pulse shapes for C{t) shown in Figure [9] we find that 
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(48) accurately predicts the total (integrated) pulse area required to make the transfer 
a-yb, i.e. /C(^)dr ^ (tr. 

An important point to note is that the time-modulated asymmetry technique 
can only be used to directly connect states with opposite parity, as is evident from 
expression (48). Like for any one-dimensional Helmholtz equation, the eigenstates of 
the Mathieu equation alternate in parity as one goes up in energy. Starting from the 
ground state, which is even, we can therefore only directly excite odd parity states. 
Furthermore, rapidly becomes small when the eigenstates a and b are far apart 
in energy, i.e. their labels j differ significantly. Crudely speaking, this is because in 
fc-space the eigenstates below the separatrix resemble those of the simple harmonic 
oscillator (Hermite polynomials) in coordinate space and are strongly peaked at the 
classical turning points The overlap integral between the initial and final state 



that occurs in (48) therefore rapidly decreases in magnitude when the classical turning 
points differ i.e. for eigenstates a and b which are far apart in energy. For this reason 
we find the surprising result that it can be far quicker to excite up to the final level by 
a series of steps via intermediate levels rather than to directly excite the upper level 
(this also allows us to excite a final state that has the same parity as the initial state). 
In Figure [9] we show the results of the stepping process 1— >2— *3^4— *5, so that 
the initial and final states have the same parity, for the case A = 10 (the separatrix 
occurs between states 5 and 6). Figure [9] also shows the pulse envelopes for each step: 
each envelope has the form a x (1 + tanh[r/6]), where a and b are constants. The total 
pulse "area" was C^R = 9-02 which is the sum of the tr values predicted by Equation 
(48), and to make the process adiabatic we set the total excitation time to be r = 400. 
In order to compare the single step with the multi step method, consider the case 
1^2^3^4^5^6, i.e. to the state immediately above the separatrix, which has 
odd parity. Equation (48) gives a total pulse "area" of 11 for the multi step method 
whereas it gives a pulse "area" of 777 for the single step 1^6. The efficiency of the 
multi step method in comparison to a single step method becomes even more striking 
as the number of steps increases, albeit at the cost of greater experimental complexity. 



7. Bragg resonances in the Josephson Junction 

In this and the final Section we consider two types of dynamics, which we shall refer to 
as Bragg scattering and Bloch oscillations, which are purely quantum effects (i.e. beyond 
mean- field) and so are not present at all in the Josephson equations ([T]) and (|2|. Rather, 
these types of motion only appear in the quantum treatment embodied by equations 
([T]) and As their names suggest, these two phenomena can be viewed as analogues 
of well-known wave scattering effects in lattices. Indeed, Haroutyunyan and Nienhuis 
[33] have previously given an analysis of the mathematical connections between atomic 
double-well systems and atomic diffraction, including the Bragg scattering analogy. Our 
purpose here is rather to present an intuitive physical discussion and refer the interested 
reader to [33] for more formal details. 
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Consider a beam of non-interacting quantum particles that form a plane wave 
incident upon a lattice whose interaction with the particles is given by the periodic 



potential V{x) = —VqcosKx, as depicted in Figure 10 Such an interaction is realized 
in experiments where atoms are diffracted from a standing wave of laser light pll |35| [36] , 
but the same basic form of refractive index also occurs in, for example, the description 
of the diffraction of light by ultrasonic waves, the original setting of the RN equations 
123]. If the total momentum of each particle is hC,, and the x-component is kg, 
the atom beam travels at angle sin^ = q/^ to the z axis and is described by the wave 
function \E'(x, z) = exp[i{^/^'^ — q^z + gx)]. The periodicity of the lattice means that it 
can only transfer momentum to the particle wave function in discrete units of h,K and 
so the diffraction of the particles is captured by the wave function [271 El] 



= exp[iy^^2 _ ^2^] ^ A„(z)exp[i(ni^ + (49) 

where Ani^z^ is the amplitude of the nth diffracted beam which travels at an angle 
tan^^ = {nK + - q^ to the z axis. The beam amplitudes are functions 

of the depth z through the diffracting medium: they describe dynamical diffraction 



in a thick grating. The wave function (49) is only an approximate description of the 



true experimental situation because it assumes that the z-component of momentum, 
namely 'h\/^ — g^, is a constant of the motion unaffected by the entry and exit from 
the periodic potential. This holds when the particles' incident energy is much greater 



than the lattice potential Yi'^ jlm ^ Vq. Substitution of (49) into the Schrodinger 
equation yields the equations for dynamical diffraction. In the paraxial approximation, 
which is valid when the diffraction angles are small enough that the term d'^An/dz'^ can 
be neglected, these equations reduce to the RN equations. The RN equations for oblique 



incidence have exactly the same form as Equation (47) that we used above to describe 
the asymmetric double- well Josephson junction, and this connection forms the basis of 
the analogy between the two cases. However, in the diffraction problem the parameters 
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that appear in (47) have a different physical origin: 

2mVo , , 

A (50) 

Ai? = 2-| (51) 

Thus, for the diffraction problem the parameter A that determines whether the system 
is in the quantum (A small) or classical (A large) regime depends upon whether the 
lattice is shallow or deep, respectively, in comparison to the recoil energy of the particles 
h'^K^/2m. The magnitude of the inbalance /S.E between wells in the double-well case 
is determined in the diffraction case by the initial transverse wavenumber q (i.e. the 
component of the initial momentum along the x-direction), or in other words, the 
angle of incidence. Finally, the dimensionless time parameter r corresponds to the 
distance z travelled through the lattice potential. In atomic diffraction experiments the 
parameters (50)-(52) can all be tuned over large ranges, something which is much harder 
to do in conventional solid state X-ray or electron diffraction experiments. For instance, 
the dimensionless Planck's constant parameter A can be controlled by the intensity of 
the laser beams forming the standing wave. However, one should not confuse in this 
discussion the interacting atoms in the double-well problem with the non-interacting 
atoms in the atomic diffraction problem since they play quite different roles in the 
models described here. 

Given the above the general treatment of the diffraction problem, we now 
specialize to Bragg scattering. Bragg resonances occur when the conditions are met 
for constructive interference of waves reflected from the different planes of a lattice and 
are therefore a pure wave phenomenon that has no counterpart for particles. Bragg 
diffraction has been observed in a number of atomic diffraction experiments e.g. |35] 
and [36]. The well known Bragg condition states that the Bragg angles 9-q satisfy 

2dsm 6b = nX (53) 

where d = 2ti / K is the period of the lattice, A = 27r/^ is the wavelength of the 
incident wave and n is a positive or negative integer. The Bragg condition (53) can 
be written in terms of the initial transverse wavenumber as ge = nK/2. The Bragg 
scattering resonance couples the incident beam (the term in Equation (49)) to the 
—nth diffracted beam (the A_n term in Equation (49)) and this corresponds to specular 
reflection from the planes of the potential. 

When the periodic potential is weak , i.e. when A is small, and we are close to a 
Bragg resonance, we can make the well known two-beam approximation and restrict our 
attention to only the Aq and A_n beams, the amplitudes of the other diffracted beams 
being much smaller PH 15^ . Taking, for example, the case where = 1, the two- 

beam solution to the RN equation (47) gives the following expressions for the intensities 
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Figure 11. A = 1. Intensity of Bragg diffracted beam as given by the analytic two-beam 

approximation \55\ , plotted as a function of AE = 2q/K (which is related to the angle of incidence via 
sin 9 = q/l;) and time r. For short times the resonance is very wide, much wider than A. For longer 
times |j4_ip varies rapidly as a function of AE and the central resonance narrows. The periodic 
oscillations in r are known as pendcUosung oscillations in X-ray diffraction theory. 



of the two beams as a function of time (equivalent to propagation depth through lattice) 



l^o(r)P 



cos 



{AE - 1)2 + A2 r/2 



(54) 



+ 



{AE - 1)2 H 
A2 



A2 



sm 



(AE - 1)2 + A2 r/2 



sm 



(AE - 1)2 + A2 r/2 



(55) 



(AE - 1)2 + A2 

We see that the width of the resonance as AE is varied is controlled by A. The 
acceptance angle (width of resonance) is wide for short times but narrows at longer 
times. This can be understood on the grounds of an energy-time uncertainty argument 
AS At > h. Notice also that the relative population of the two beams oscillates as a 



function of the time/depth through the grating r with a period 27r/y (AE — 1)2 + A2. 
In the field of electron diffraction these oscillations are known as "pendellosung" 



(pendulum solutions). The cases where 



\n\ 



2,3,4... are more complicated but 



analytic expressions can be obtained by adiabatically eliminating the intermediate 
amplitudes so that one still has a two-beam solution [37] . 

The analytic result in the two-beam approximation (55) for the intensity of the first 

P is plotted in Figure 11 Meanwhile, Figure 12 displays the 



Bragg diffracted beam \A. 
results of an exact numerical solution of the full set of time-dependent RN equations (47). 
As expected, we find that the approximate analytical results agree with the numerical 
ones providing A is small and we are close to a Bragg angle. In particular, the top row of 
graphs in Figure 12 show the temporal evolution of the two strongly-coupled beams at 
the Bragg resonances (a) AE = 1 and (b) AE = 2. The pendellosung are clearly visible, 
particularly in (a). In the bottom row we see the classic Bragg resonance structure as a 
function of the incident angle AE. The value of r chosen for the bottom row of Figure 



12 is such that the population of the relevant Bragg scattered beam is at a peak, i.e. 



half a pendellosung period. 
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Figure 12. A = 1. Upper panels: Pendellosung oscillations found by numerical solution of 
the full RN equations | |47| . In the diffraction problem the pendellosung are between two diffracted 
beams as a function of depth through the lattice. In the Josephson problem the pendellosung 
become oscillations between two number difference states as a function of time. Which two diffracted 
beams (number states) are involved depends on the value of the incident angle (energy asymmetry) 
AE = 2q/K. In (a) AE = 1 and in (b) AE = 2. All the other beams (states) which are not 
shown have much smaller populations. Lower panels: the coupling between the two number difference 
eigenstates (diffracted beams) has a clear resonance structure as a function of AE. Resonances occur 
at AE = . . . , —2, —1,1,2,.... The lower panels are shown at times approximately equal to half their 
respective pendellosung periods: c) r = 3.4045 (d) r = 7.1525. 



We can now use these simple results from the theory of the diffraction of non- 
interacting particles to predict phenomena for the case of interacting particles in a 
double- well potential. Recalling that, according to Equation ([9]), the amplitudes Aq and 
A_n refer to number difference states we see that the analogue of a Bragg resonance 
in the Josephson problem is a tunnelling resonance between two number difference, i.e. 
the swapping of a precise number of atoms between the two wells [33] as the amplitudes 
Aq and A_n oscillate. In order to achieve the analogous initial conditions as in the 
"real" Bragg scattering scenario described above, where the incident atom beam was 
well coUimated so that at r = we have Aq = 1 and A„^o = 0, it is necessary to start in 
a single population difference state. For simplicity, in what follows we shall assume that 
the initial state is the precisely balanced state with k = (this is the most likely result in 
a balanced double- well) and so the initial conditions are Aq = 1 and An^Q = 0. We shall 
return at the end of the Section [8] to consider the challenging demands this places on an 
experimental realization. At time r = we suddenly switch on an energy asymmetry 
AE which is held at a constant value. The temporal evolution of the amplitudes is 



then exactly that shown in Figures 11 and 12 As shown in the top row of Figure [12 



Classical versus quantum dynamics of the atomic Josephson junction 



23 



by tuning the energy asymmetry AE to a resonance we find pendellosung between two 
number difference states. At resonance, the Josephson pendellosung oscillations between 
the states n = and n = —1 occur with a period 

'^pendellosung = 27r/A . (56) 
This has a different dependence upon A than the Josephson oscillation period (as given 



by Equation (37)) and, combined with the resonance behaviour, would provide an 
experimental signature of the Bragg analogue. Although we have only plotted the 
\n\ = 1 and \n\ = 2 Bragg resonances, corresponding, respectively, to the tunnelling 
oscillation of one and two atoms between the wells, it is possible to isolate resonances 
involving the precise transfer of larger numbers of atoms. 

8. Band Structure, Bloch oscillations and Schrodinger cats in the 
Josephson Junction 

It is instructive to re-cast our discussion of diffraction in a lattice in terms of band 
theory [55]. We shall see that this gives a deeper understanding of the Bragg resonance 
phenomena in a double-well and also leads naturally to analogies to a class of adiabatic 
phenomena related to Bloch oscillations. According to the Bloch theorem, the stationary 
eigenfunctions for a quantum particle in the lattice potential V{x) = —VqcosKx can 
be factorized 

?/^«'^(x) = exp(iga;)f/^'^(2;) (57) 

in terms of a plane wave part exp(igx) and a spatially dependent amplitude part Lf'^i^x) 
which is periodic with the same period as the lattice: U'^'^{x + 2tt/K) = W'^lx). 
The wave function ip'^'^{x) depends on two quantum numbers: q which gives the 
quasimomentum hq, and the band index j. In the diffraction problem presented above 
the quasimomentum is nothing but the x-component of the momentum of the incident 
wave. The band index j is also already familiar to us: we use the same symbol we used 
to label the eigenenergies and eigenvectors of the Josephson problem because it has an 
identical meaning here. However, because in the Josephson problem the wave function 
must be 27r-periodic, there we had to set q = 0, thus seemingly loosing the full richness of 
the general problem of a quantum particle in a periodic potential. Substituting ip^'^{x) 
into the time-independent Schrodinger equation yields 

f/'?'^ (x) - Vo cosKx U'^'^ix) = e'^'^U'^^^ix). (58) 




The eigenenergy e'^'\ which is a function of both the band index j and the 



quasimomentum q, has the well known band structure shown in Figure 13 The 
periodicity of U'^'^{x) means that we can expand it as a Fourier series 

U'''^^) = ^^ E Bt:^ e^^{;mKx). (59) 
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(a) (b) 

Figure 13. Energy bands of a quantum particle in a sinusoidal potential as a function of the 
dimcnsionlcss quasimomentum AE = 2q/K. (a) Shows the energies of the relevant momentum basis 
states in the absence of the sinusoidal potential and (b) shows the actual energy bands in the presence 
of the sinusoidal potential with A = 0.2. In particular, the thick red parabola in (a) is the dispersion 
relation of a free particle. This free particle corresponds to the incident wave in the diffraction 
problem discussed in Section [7] The other parabolas also represent free particles but their momenta 
have been shifted by an integer number of units of hK. Together these states form the natural basis 
states inside the sinusoidal potential (see text). The energy level structure shown in (b) is calculated 
using the Raman-Nath equations ( |60| l. We see that the parabolic basis states are coupled such that 
their intersections become avoided crossings. Bragg resonances lie at the places where the thick red 
parabola crosses the other parabolas. Bloch oscillations correspond to an adiabatic evolution along a 
single band as AE is varied. 



Substituting this into (58) yields 



n 



B: 



A 
2 



1,3 



)B 



(60) 



In this formula A is the parameter defined in Equation (50), and the dimensionless 

2me'^'^ / Ti^ K'^ . Comparing (60) with the RN Equations 

(a) 



energy £'^'^ is given by E"^'- 



(29) for the asymmetric double- well, we see that the two are identical if: (a) we assign 
2q/K = AE, just as we already did in (51), i.e. we recognize that the tilt AE is really 
equivalent to a quasimomentum, and (b) we put 

f_ . AE^ 



gq,3 



Ei 



(61) 



Recall that E^^ is the eigenvalue corresponding to the jth eigenvector of the asymmetric 
double- well and so is a function of AE. 

The fact that the tilt asymmetry AE between the wells allows us to introduce a 
term that plays the role of a quasimomentum into the Josephson problem even though 
the wave function is strictly 27r-periodic is of considerable significance; it is crucial 
for many of the phenomena discussed so far and also for Bragg resonances and Bloch 
oscillations (see below). In the Josephson problem the quasimomentum is introduced 
directly into the hamiltonian (see Equation (23)) rather than via the wave function: it 
is a parameter rather than a dynamical variable. 

An example of the band structure of the eigenenergies £'^'^ obtained by solving the 
RN equations (60 ) is shown in Figure 13 In the absence of the periodic potential (A = 0) 
the eigenstates are exp(iga;) corresponding to a free particle. Their dispersion relation 
h'^q^/2m is plotted as the thick red parabola in Figure 13 a). Also shown in Figure 



Classical versus quantum dynamics of the atomic Josephson junction 



25 



13 a) are the displaced parabolas fi^iq + nKy/2m centered at the positions q = —uK 



which are the dispersion relations of the momentum eigenstates exp[i(g + nK)x\, where 
n = 0, ±1, ±2, ±3 . . .. These are the basis states used in the expansion (59) of the wave 
function inside in the periodic potential. Inside the periodic potential the momentum 
basis states are coupled, which is the physical content of the the RN equation (60). For 
shallow lattices (A < 1) this coupling is weak and is only significant close to the points 
where the free parabolas would otherwise have crossed and caused a degeneracy. The 
effect of the coupling is to turn the degeneracies into avoided crossings, thereby forming 
continuous energy bands as a function of g, separated by gaps of forbidden energies, 
as can be seen in Figure 13 b). Technically, each band is a superposition of an infinite 
number of basis states although when the potential is weak only a small number make a 
significant contribution at any particular value of q. However, the particular combination 
required depends on q. Take, for example, the first Brillouin zone (BZ) defined as the 
region —K/2 < q < K/2. When A is small the first band in the first BZ can be 
adequately described using only the n = state exp[igx] and the two |n| = 1 states 
exp[i(g — K)x\ and exp[i(g + K)x\, which are needed to take into account the coupling 
between the free particle states at the edges of the first BZ. For instance, at the right 
hand edge of the first BZ, where q = K/2, there is a strong coupling between between 
the n = and n = —1 free particle states exp[iga:] and exp[i(g — K)x], respectively, 
whose energies would otherwise be degenerate. Using degenerate perturbation theory 
one finds that the band splitting at the avoided crossing between the first and second 
bands is given by 

Ei=' - = A. (62) 

On the left hand side of the first Brillouin zone at g = —K/2 it is the n = and n = 1 
free particle states exp[igx] and exp[i(g + K)x], respectively, that are strongly coupled. 

Bragg resonances have a simple interpretation in this energy band picture: they 
correspond to conservation of energy and momentum. Referring to Figure [l3|(a), 
Bragg scattering takes place where the dispersion relation of the incident wave 
(thick red parabola) crosses those of the other basis states. We see that Bragg 
resonances correspond to having an incident quasimomentum given by ge = {n/2)K 
or, equivalently, AE^ = n where n is a postive or negative integer. Bragg scattering 
therefore takes place at the edges of the Brillouin zones. Consider a wave exp[iga;] 
(ignoring the trivial dependence on z) incident somewhere in the first BZ. If it enters 
with a value of q not close to g = ±i^/2 then, according to the above discussion, when 
A is small this wave coincides with the eigenfunction describing the first band and it 
propagates through the lattice unmodified and hence undefiected. The weak potential 
is not capable of coupling the wave to higher bands unless the wave enters at a Bragg 
angle. In the case that the wave is incident close to either of the first Bragg angles 
AE = ±1 the wave function exp[igx] of the incident wave is seen to be a superposition 
of the first and second band eigenfunctions. The beating between the two bands leads 
to the pendellosung. Similarly, at AE = ±2 there are Bragg resonances between the 



Classical versus quantum dynamics of the atomic Josephson junction 



26 




-4 -3 -2 -1 1 2 3 4 -4 -3 -2 -1 1 2 3 4 

AE AE 



(a) (b) 

Figure 14. Energy band structure of the Josephson double-well system. In the Josephson double- 
well problem the role of the quasimomentum is played by the tilt asymmetry AE between the double- 
wells. Each line in (a) corresponds to the energy of a particular number difference eigenstate as given 
by Equation ( |63| . These states are the eigensolutions of the double- well problem in the absence of 
tunnelling (A = 0). (b) shows the actual energy level structure in the presence of tunnelling (A = 0.2). 
Note that every apparent level crossing in (b) is actually an avoided crossing. The bottom curve is 
the ground band i?^, the next up is for E^, the next is for E^ etc. The eigenenergies Ei of the 
asymmetric double-well for a particular value of AE correspond to vertical slices through this band 
structure. The thick red horizontal line in (a) is the energy of the ra = number difference eigenstate 
and is equivalent to the thick red parabola shown in Figure [Tsja). In fact, using the transformation 
between f and we find that the entire band structure shown here is equivalent to that 
shown in Figure [13] 



incident wave and the second and third bands and so on for larger angles of incidence. 
The magnitude of the band gap at the avoided crossings gives an indication of the 
strength of the corresponding Bragg resonance. Generalizing (62) we find that the band 
gap between the j + 1 and the jth bands scales as A-' and so the higher resonances 
become weaker when A < 1. In the case of atomic diffraction this can be understood 
physically by noting that 2j photons must be exchanged between the atoms and the 
laser beams at the jth Bragg resonance. 

The foregoing analysis suggests a band structure interpretation of tunnelling 
resonances in the Josephson problem. The energy eigenvalues of the RN equations 
(29) as a function of the tilt asymmetry AE between the wells have the structure shown 
in Figure 14 In the quantum regime A is small meaning that the interaction energy 
Ec dominates the tunnelling energy Ej. In this limit it makes sense to choose the 
basis states to be the number difference eigenstates exp[in(/)], which are analogous to the 
momentum eigenstates we used in the diffraction problem, and indeed the RN equations 
are written in this basis. However, unlike in the diffraction problem, the quasimomentum 
phase factors exp[igx] = exp[iA£'/2 0] are not included in the basis states because they 
occur in the hamiltonian not the wave function, which must be 27r-periodic as noted 
above. Furthermore, the energies associated with the number difference eigenstates are 
not parabolic but are linear in AE 

E{n,AE)=n^ + nAE. (63) 

These energies are shown as the straight lines in Figure [l4|( a). As can be seen in Figure 
[T4Fb), there is a band structure associated with the Josephson double- well problem, but 
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it is distorted by gravity. When we plot the bands in terms of their "true" energies 
El (rather than i^^'-') the forbidden energy gaps in the spectrum disappear except for 
one, just below E = 0, although this itself becomes extremely narrow for large AE, and 
in any case our hamiltonian eventually breaks down at large tilts where the condition 
AE -C is violated. However, with the trivial redefinition (61) of the energy from Ei 
to S'^'^ the standard band structure shown in Figure 13 is recovered. 

Let us examine in more detail the band structure of the Josephson double-well 
problem as illustrated in Figure 14 for small A. For small tilts {AE ^ 1) the first band 
is adequately described by the number difference eigenfunction exp[iO0] corresponding 
to zero atom number difference between the wells and represented by the thick red 
horizontal line in Figure 14 a). Similarly, the second band in this range consists of a 
superposition of the states exp[i0] and exp[— 10], which describe cases where one particle 
has been transferred to the left side and one has been transferred to the right side, 
respectively. Near the Bragg resonance at the right hand edge of the first BZ where 
AE = 1, the states exp[iO0] and exp[— 10] are strongly coupled leading to the avoided 
crossing shown in Figure [T4|^b). By analogy with the diffraction case, the first Josephson 
Bragg resonance occurs when the double-well system is suddenly tilted from AE = to 
AE = 1. This projects the initial n = state over the superpositions of the n = and 
n = 1 states that make up the eigenfunctions giving the first two bands at AE = 1. A 
single atom will then oscillate back and forth between the two wells. 

The successively higher Bragg resonances, i.e. the places where the thick red line 
in Figure [l4|a) crosses the lines located successively further away from the origin, 
correspond to places where the n = balanced state couples strongly to states with 
successively higher values of n. When the high n state is macroscopically large, n = P, 
say, where P ^ 1, the system oscillates between having a population difference of 
precisely zero and precisely P atoms. Denoting by \n) the number difference kets, 
at a Bragg resonance the double-well is in a superposition of the symmetric and 
antisymmetric Schrodinger cat states \E'|^ = (|0) ± |P))/-\/2, which are themselves 
superpositions of macroscopically distinguishable states. The lower of the two bands 
making up the avoided crossing corresponds to the symmetric state and the higher to 
the antisymmetric state. In order to realize a single Schrodinger cat state, as opposed 
to a superposition, i.e. either or ^f^*", it is necessary to occupy just one band. One 
way to achieve this is via Bloch oscillations rather than Bragg resonances. 

Bloch oscillations are in a sense the conjugate phenomenon to Bragg scattering 
because they rely on adiabatic evolution rather than sudden projection. Conventionally, 
Bloch oscillations occur when a constant external force F is applied to a quantum 
particle in a periodic potential, e.g. an electron in a crystal lattice subject to a constant 
electric field. It transpires that in a periodic potential it is the quasimomentum which 
obeys Newton's second law under the action of the force 

qit) = go - Ft/h (64) 
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a result known as Bloch's acceleration theorem |39lll0lllT]. In this formula go is the value 
of g at t = 0. Thus, under the influence of the external force the quasimomentum travels 
along the band at a linear rate, and if this motion is adiabatic the system remains in 
the same band. The periodicity of the band structure means that rather than uniformly 
accelerating the particle undergoes an oscillatory motion. Non-adiabatic corrections to 
this evolution involve Landau-Zener transitions to other bands at the avoided crossings. 
By contrast, Bragg scattering is inherently non-adiabatic and involves suddenly exciting 
a superposition of two bands. 

In the case of the double-well Josephson problem an analogy to Bloch oscillations 
occurs when the tilt asymmetry is slowly changed. Imagine starting at AE = in the 
first band in Figure 14 and slowly increasing AE, thereby sweeping the system along 
the first band. The initial state therefore has the same number of atoms in each well. 
If the sweep proceeds adiabatically the system remains in the first band. At the first 
avoided crossing it smoothly evolves into the new ground state so that it is now in a 
state where exactly one particle has been transferred into the lower well. If the adiabatic 
sweep continues we move successively through states with more and more atoms in the 
lower well. This can continue until all the particles are in the lower well. In Figure 
[l4|^b) it appears that the avoided crossings between the first and second band become 
smaller and smaller further from the origin, and so the sweep must become slower and 
slower to remain adiabatic, but this is an illusion. From the transformation (61) we see 
that the bands gaps are identical to those in Figure [l3|(b) where one can see that they 
are independent of AE. Note, however, that our approximate hamiltonian (23) is only 
valid as long as AE <^ N. When this condition is violated the band gaps will become 
a function of AE and so the adiabaticity condition will in general depend on AE for 
large asymmetries. 

If the avoided crossing traversals are not entirely adiabatic we mix in some Bragg 
scattering-like character into the evolution. Consider a non-adiabatic traversal of the 
avoided crossing between the first BZ and second BZ: the system is then put into a 
superposition of the first two bands (i.e. a superposition of the balanced state and the 
state with one particle transferred) and the system will be set into oscillation. Non- 
adiabatic traversals imply that the final state when AE becomes large will not be the 
ground state with all the atoms in one well but rather an excited state with atoms 
oscillating between the wells. 

If one wishes to excite a single Schrodinger cat state \E'^*" in the double-well it is 
necessary to use the time oscillating asymmetry method of exciting a single eigenstate 
as described in Section |6j Consider starting, as before, in the precisely balanced state 
|0) at AE = 0. If one were to simply increase AE adiabatically then at the first avoided 
crossing one generates the state (|0) + | — l))/-\/2- Remaining in this first band only 
generates superpositions between the successive neighbouring parabolas shown in Figure 
[Tsj^a) which only differ by a particle number difference of a single particle. To generate 
a superposition of two states differing by a macroscopically large particle number one 
should use the time oscillating asymmetry to excite a higher band. In fact, remaining 
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on the AE = axis generates superpositions of the form {\n) + | — n))/\^. However, if 
one tries to directly excite to an avoided crossing where two bands nearly touch there 
is a considerable risk of exciting both bands. Therefore, a superior method would be 
to first adiabatically sweep the tilt asymmetry whilst still in the first band to the point 
AE = 0.5 where the two higher lying bands are maximally separated, apply the time 
oscillating tilt to excite just one of them, then adiabatically sweep AE back to zero to 
reach the high lying avoided crossing. 

Finally, we briefly discuss some points concerning the realization of Josephson Bragg 
scattering and adiabatic phenomena such as Bloch oscillations. Bragg scattering requires 
that the initial state is a single population difference state. One way this can be achieved 
is if the double well begins in its ground state in the quantum limit where A <^ 1 (low 
tunnelling rate regime). Then, if the tilt is suddenly changed to one of the Bragg angles 
AEb = n the initial state is projected over the two bands with which it is resonant (see 



Figure 13). Because the Bragg resonances have a finite width there is some tolerance to 
errors in AE. Higher Bragg resonances involving the tunnelling of a large but precise 
number of atoms are harder to achieve because the higher resonances are weaker when 
A < 1 as mentioned above. Going to larger A has the effect of mixing in a larger number 
of population difference basis states into each band. This means that one can no longer 
claim, for example, that the ground band at AE = solely consists of the n = state. 
When it comes to Bloch oscillations driven by a sweep in AE, adiabaticity is most likely 
to be maintained in the semiclassical (meanfield) limit of large A where the tunnelling 
ensures a large splitting between states at the avoided crossings. However, if one is 
interested in achieving a precise number difference of atoms between the two wells using 
this adiabatic method then, for the reasons already mentioned, it pays to have a small 
value of A. All of these considerations need to be set in the context of the very significant 
experimental challenge of cooling the double well system to its ground state [21]. This 
is the desirable initial state for demonstrating Bragg scattering, Bloch oscillations and 
also Schrodinger cat states. The robustness of these phenomena to finite temperature 
effects will be the subject of future work. 

9. Conclusions 

We have analysed the atomic Josephson junction from the point of view that the 



parameter 1/vA = y Ec/{2Ej) is proportional to Planck's constant. This provides 
a simple way to predict when the dynamics obeys the classical (mean-field) Gross- 
Pitaevskii theory and when it must be quantised. If the system is set into motion 
by suddenly removing an asymmetry between the wells both low energy Josephson 
plasmons and high energy rotor excitations are accessible, depending on the magnitude 
of the asymmetry. For small A the system is very quantum and both excitations rapidly 
deviate from the classical mean-field predictions, the former by a change in frequency 
and the latter by undergoing collapses and revivals. We give an expression for the period 



of the revivals, see Equation (45) and the surrounding discussion. 
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Between the low energy and high energy regimes hes a classical separatrix. We 
find that the quantum and classical evolution always quickly diverge from each other 
for motion close to the separatrix. This is true even when in other respects the system 
is expected to behave classically, i.e. A is large. We interpret this divergence as being 
due to quantum tunnelling and show that quantum and classical must deviate after one 
quarter of the classical period. 

A second method for exciting the system is to have a tilt asymmetry which is 
periodically modulated in time. By tuning the modulation frequency this method allows 
the excitation of a single eigenstate which is of course a non-classical state. In the 
quantum limit A — > this method provides a way of generating states of the form 
(|n) ± I — n))/\/2 where n is the difference in the number of particles between the two 
wells. When n becomes large these are Schrodinger cat states involving macroscopically 
distinguishable superpositions. 

Finally, we have discussed at length the analogy between the asymmetric (tilted) 
double-well problem and the diffraction of waves by a periodic lattice, including an 
analysis in terms of band structure familiar from solid state physics. Bragg scattering 
in the diffraction problem corresponds to tunnelling resonances in the Josephson problem 
where a precise number of atoms oscillate between the wells in close analogy with the 
pendellosung oscillations. These resonances are not present in the classical Gross- 
Pitaevskii theory. The role of the quasimomcntum in the band structure analysis is 
played by the energy difference between the two wells. Bragg resonances correspond to 
a sudden tilt of the double wells to a precise final value (the "Bragg angle" ) whereas a 
slow steady increase of the tilt asymmetry is analogous to Bloch oscillations. 
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